library(readxl)     
library(dplyr)      
library(lubridate)  
library(writexl)    

ts <- sprintf("2025-%02d", 10:12)

load("~/Salidas/c. v_02.RData")
load("~/Salidas/c. Estimout_01_06.RData")
load("~/Salidas/b. MDATA4.RData")
DATA <- DATA %>% filter(EMPRESA != "VIRGIN")
EMPRESA <- DATA$EMPRESA

i <- which(DATA$t %in% ts)
v <- v[i, ]
delta <- delta[i]
cm <- cm[i]
EMPRESA <- EMPRESA[i]

DATA <- DATA %>%
  slice(i) %>%
  mutate(
    delta = delta,
    cm = cm
  )

p <- DATA$p
S <- DATA$s
t <- DATA$t
NS <- ncol(v)
mkti <- match(DATA$t, unique(DATA$t))

DATA$deltar <- DATA$delta + zeta[2] * p 

dx_obs <- DATA %>%
  mutate(firm_id = EMPRESA) 

H_obs <- matrix(0, nrow = nrow(DATA), ncol = nrow(DATA))

for(m_idx in seq_along(unique(t))) {
  market_rows <- which(mkti == m_idx)
  firms_in_market <- dx_obs$firm_id[market_rows]
  matching_firms <- outer(firms_in_market, firms_in_market, "==") * 1
  H_obs[market_rows, market_rows] <- matching_firms
}

p0 <- p
pup <- 0.5
O <- nlminb(start = p, 
            objective = function(x) sum(mktsimul(x, p0 = p0, 
                    deltar = DATA$deltar, 
                    cm = DATA$cm, 
                    H = H_obs, 
                    v = v, 
                    zeta = zeta, 
                    pup = 0.5)^2),
            lower = p0 * (1 - pup),
            upper = p0 * (1 + pup))
p_simul <- O$par
s_simul <- share_2(p = p_simul, deltar = DATA$deltar, v = v, zeta = zeta) 

DATA <- DATA %>%
  mutate(
    EMPRESA = recode(EMPRESA, "MOVISTAR" = "TIGO") 
  ) %>%
  group_by(t, EMPRESA, posp) %>%
  mutate(
    cm = if_else(EMPRESA == "TIGO", 
                 0.8814 * weighted.mean(cm, s, na.rm = TRUE), 
                 cm), 
    
    deltar = if_else(EMPRESA == "TIGO", 
                 weighted.mean(deltar, s, na.rm = TRUE), 
                    deltar) 
  ) %>% 
  ungroup()

dx <- DATA %>%
  mutate(firm_id = ifelse(EMPRESA %in% c("TIGO", "MOVISTAR"), "TIGO", EMPRESA))

H_Int <- matrix(0, nrow = nrow(DATA), ncol = nrow(DATA))

for(m_idx in seq_along(unique(t))) {
  market_rows <- which(mkti == m_idx)
  firms_in_market <- dx$firm_id[market_rows]
  matching_firms <- outer(firms_in_market, firms_in_market, "==") * 1
  H_Int[market_rows, market_rows] <- matching_firms
}

p0 <- p
pup <- 0.5
I <- nlminb(start = p, 
            objective = function(x) sum(mktsimul(x, p0 = p0, 
                    deltar = DATA$deltar, 
                    cm = DATA$cm, 
                    H = H_Int, 
                    v = v, 
                    zeta = zeta, 
                    pup = 0.5)^2),
            lower = p0 * (1 - pup),
            upper = p0 * (1 + pup))
P_Int <- I$par
S_Int <- share_2(p = P_Int, deltar = DATA$deltar, v = v, zeta = zeta) 

dx <- DATA %>%
  mutate(firm_id = ifelse(EMPRESA %in% c("COMCEL", "TIGO", "MOVISTAR"), "CARTEL", EMPRESA))

H_coor <- matrix(0, nrow = nrow(DATA), ncol = nrow(DATA))

for(m_idx in seq_along(unique(t))) {
  market_rows <- which(mkti == m_idx)
  firms_in_market <- dx$firm_id[market_rows]
  matching_firms <- outer(firms_in_market, firms_in_market, "==") * 1
  H_coor[market_rows, market_rows] <- matching_firms
}

p0 <- p
pup <- 0.5
C <- nlminb(start = P_Int, 
            objective = function(x) sum(mktsimul(x, p0 = p0, 
                    deltar = DATA$deltar, 
                    cm = DATA$cm, 
                    H = H_coor, 
                    v = v, 
                    zeta = zeta, 
                    pup = 0.5)^2),
            lower = p0 * (1 - pup),
            upper = p0 * (1 + pup))
P_coor <- C$par
S_coor <- share_2(p = P_coor, deltar = DATA$deltar, v = v, zeta = zeta) 

a <- data.frame(
  FECHA           = DATA$t, 
  EMPRESA         = DATA$EMPRESA,
  POST            = DATA$posp,
  deltar          = DATA$deltar,
  M               = DATA$M, 
  CM              = DATA$cm,
  P               = DATA$p,
  P_SIMULADO      = p_simul,
  P_INTEGRACION   = P_Int, 
  P_COORDINACION  = P_coor,
  S               = DATA$s,
  S_SIMULADO      = s_simul, 
  S_INTEGRACION   = S_Int, 
  S_COORDINACION  = S_coor,
  stringsAsFactors = FALSE
)

firm_vector <- a$EMPRESA
target_firm <- "TIGO"
idx_change <- which(firm_vector == target_firm)
idx_fixed  <- which(firm_vector != target_firm)
p_optim0 <- a$P_COORDINACION[idx_change]
p_fixed <- a$P_COORDINACION[idx_fixed]

iter_count <<- 0
T_1 <- optim(par = p_optim0, 
             fn = best_response_2, 
             method = "L-BFGS-B",
             control = list(fnscale = -1, trace = 1))

a$P_TRAICION <- numeric(nrow(a))
a$P_TRAICION[idx_change] <- T_1$par
a$P_TRAICION[idx_fixed]  <- p_fixed

mui    <- - zeta[1]*v*a$P_TRAICION
deltar <- a$deltar
delta  <- deltar - zeta[2]*a$P_TRAICION
S_TRAICION   <- share(delta = delta, mui = mui)
a$S_TRAICION <- S_TRAICION$s

firm_vector <- a$EMPRESA
target_firm <- "COMCEL"
idx_change <- which(firm_vector == target_firm)
idx_fixed  <- which(firm_vector != target_firm)
p_optim0 <- a$P_COORDINACION[idx_change]
p_fixed <- a$P_COORDINACION[idx_fixed]

iter_count <<- 0
T_2 <- optim(par = p_optim0, 
             fn = best_response_2, 
             method = "L-BFGS-B",
             control = list(fnscale = -1, trace = 1))

a$P_TRAICION_2 <- numeric(nrow(a))
a$P_TRAICION_2[idx_change] <- T_2$par
a$P_TRAICION_2[idx_fixed]  <- p_fixed

mui    <- - zeta[1]*v*a$P_TRAICION_2
deltar <- a$deltar
delta  <- deltar - zeta[2]*a$P_TRAICION_2
S_TRAICION_2   <- share(delta = delta, mui = mui)
a$S_TRAICION_2 <- S_TRAICION_2$s

a <- a %>%
  group_by(FECHA, EMPRESA) %>%
  mutate(
    Pi = sum((P-CM)*S*M),
    Pi_INTEGRACION  = sum((P_INTEGRACION-CM)*S_INTEGRACION*M),
    Pi_COORDINACION = sum((P_COORDINACION-CM)*S_COORDINACION*M),
    Pi_TRAICION     = sum((P_TRAICION-CM)*S_TRAICION*M),
    Pi_TRAICION_2   = sum((P_TRAICION_2-CM)*S_TRAICION_2*M)
  ) %>%
  ungroup()

b <- a %>%
  group_by(FECHA, EMPRESA) %>%
  summarise(
    P              = weighted.mean(P, S, na.rm = TRUE),
    P_SIMULADO     = weighted.mean(P_SIMULADO, S_SIMULADO, na.rm = TRUE),
    P_INTEGRACION  = weighted.mean(P_INTEGRACION, S_INTEGRACION, na.rm = TRUE),
    P_COORDINACION = weighted.mean(P_COORDINACION, S_COORDINACION, na.rm = TRUE),
    P_TRAICION     = weighted.mean(P_TRAICION, S_TRAICION, na.rm = TRUE),
    P_TRAICION_2   = weighted.mean(P_TRAICION_2, S_TRAICION_2, na.rm = TRUE),
    S              = sum(S, na.rm = TRUE),
    S_SIMULADO     = sum(S_SIMULADO, na.rm = TRUE),
    S_INTEGRACION  = sum(S_INTEGRACION, na.rm = TRUE),
    S_COORDINACION = sum(S_COORDINACION, na.rm = TRUE),
    S_TRAICION     = sum(S_TRAICION, na.rm = TRUE),
    S_TRAICION_2   = sum(S_TRAICION_2, na.rm = TRUE),
    Pi             = mean(Pi, na.rm = TRUE),
    Pi_INTEGRACION = mean(Pi_INTEGRACION, na.rm = TRUE),
    Pi_COORDINACION= mean(Pi_COORDINACION, na.rm = TRUE),
    Pi_TRAICION    = mean(Pi_TRAICION, na.rm = TRUE),
    Pi_TRAICION_2  = mean(Pi_TRAICION_2, na.rm = TRUE),
    .groups = 'drop' 
  )

c <- b %>%
  group_by(EMPRESA) %>%
  summarise(
    P              = mean(P, na.rm = TRUE),
    P_SIMULADO     = mean(P_SIMULADO, na.rm = TRUE),
    P_INTEGRACION  = mean(P_INTEGRACION, na.rm = TRUE),
    P_COORDINACION = mean(P_COORDINACION, na.rm = TRUE),
    P_TRAICION     = mean(P_TRAICION, na.rm = TRUE),
    P_TRAICION_2   = mean(P_TRAICION_2, na.rm = TRUE),
    S              = mean(S, na.rm = TRUE),
    S_SIMULADO     = mean(S_SIMULADO, na.rm = TRUE),
    S_INTEGRACION  = mean(S_INTEGRACION, na.rm = TRUE),
    S_COORDINACION = mean(S_COORDINACION, na.rm = TRUE),
    S_TRAICION     = mean(S_TRAICION, na.rm = TRUE),
    S_TRAICION_2   = mean(S_TRAICION_2, na.rm = TRUE),
    Pi             = sum(Pi, na.rm = TRUE),
    Pi_INTEGRACION = sum(Pi_INTEGRACION, na.rm = TRUE),
    Pi_COORDINACION= sum(Pi_COORDINACION, na.rm = TRUE),
    Pi_TRAICION    = sum(Pi_TRAICION, na.rm = TRUE),
    Pi_TRAICION_2  = sum(Pi_TRAICION_2, na.rm = TRUE),
    .groups = 'drop' 
  )

a$W              <- welfare(p = a$P, deltar = a$deltar, v = v, zeta = zeta, M = unique(a$M))
a$W_INTEGRACION  <- welfare(p = a$P_INTEGRACION, deltar = a$deltar, v = v, zeta = zeta, M = unique(a$M))
a$W_COORDINACION <- welfare(p = a$P_COORDINACION, deltar = a$deltar, v = v, zeta = zeta, M = unique(a$M))

a$DELTA   <- (a$Pi_TRAICION - a$Pi_COORDINACION) / (a$Pi_TRAICION - a$Pi_INTEGRACION)
a$DELTA_2 <- (a$Pi_TRAICION_2 - a$Pi_COORDINACION) / (a$Pi_TRAICION_2 - a$Pi_INTEGRACION)

CM <- data.frame(FECHA = a$FECHA, EMPRESA = EMPRESA, CM = cm)

Resultados <- list(
  "CM" = CM,
  "DESAGREGADO" = a,
  "MENSUAL" = b,
  "TOTAL" = c
)

write_xlsx(Resultados, "c. Simulación.xlsx")
